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Abstract 

This paper aims to mathematically advance the field of quantitative thermo-acoustic 
imaging. Given several electromagnetic data sets, we establish for the first time an 
analytical formula for reconstructing the absorption coefficient from thermal energy 
measurements. Since the formula involves derivatives of the given data up to the third 
order, it is unstable in the sense that small measurement noises may cause large errors. 
However, in the presence of measurement noise, the obtained formula, together with 
a noise regularization technique, provides a good initial guess for the true absorption 
coefficient. We finally correct the errors by deriving a reconstruction formula based on 
the least square solution of an optimal control problem and prove that this optimization 
step reduces the errors occurring and enhances the resolution. 
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1 Introduction 

Hybrid imaging modalities are based on a multi-wave concept. Different physical types of 
waves are combined into one tomographic process to alleviate deficiencies of each separate 
type of waves, while combining their strengths. Multi-wave systems are capable of high- 
resolution and high-contrast imaging [1, 17]. Quantitative thermo-acoustic tomography is 
an emerging hybrid modality [14, 12]. It allows to determine the absorption distribution of 
a tissue from boundary measurements of the pressure induced by electromagnetic heating. 
Other examples of hybrid modalities are acousto-electric tomography [3, 2, 6, 9, 13, 21, 32, 
33], magnetic resonance electrical impedance tomography [20, 28, 26], magnetic resonance 
elastography [8, 25, 23], impedance-acoustic tomography [18], photo-acoustic [31, 22, 4], 
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quantitative photo-acoustic tomography [5, 11, 27], magneto-acoustic imaging [7], and vibro- 
acoustography [16]. 

The aims of this paper are to derive an exact formula for the absorption coefficient 
from noiseless thermo-acoustic measurements and to correct the errors of in the presence 
of measurement noise. The former task is motivated by the knowledge of the ratio between 
two modified data. For the latter purpose, we show how to regularize the exact formula 
and propose an optimal control algorithm to achieve a resolved image starting from the 
regularized one. As far as we know, our exact formula in this paper together with the one 
successfully derived in [6] are among a few exact formulas in hybrid imaging. Moreover, the 
fine analysis of the effect of measurement noise on the image quality and the proof that an 
optimal control approach starting from the regularized images yields a resolved one have 
never been done elsewhere. 

To describe our approach, we employ several notations. Let X be a smooth bounded 
domain in M. d , d = 2 or 3. Let dX denote the boundary of X and let v be the outward 
normal at dX. For m a non- negative integer, we define the space H m (X) as the family of all 
m times weakly differentiable functions in L 2 (X) , whose weak derivatives of orders up to m 
are functions in L 2 (X). We let H^(X) be the closure of C™(X) in H m (X), where C™(X) 
is the set of all infinitely differentiable functions with compact supports in X. Finally, we 
introduce the space i? 1 / 2 (9X) of traces on dX of all functions in H X (X). 

Let q be a positive real-valued function on X. Consider the Helmholtz problem: 

(A + k 2 + ikq)u = 0, 

(1.1) 

v ■ Vu — iku = g, x G dX, 

which is the scalar approximation of Maxwell's equations. Here, k > is the wave number, 
g is a boundary datum, and u is the electrical field. The Robin boundary condition approx- 
imates Sommerfeld's radiation condition at high frequencies [15, 19]. For simplicity, instead 
of considering the Helmholtz equation on the whole Euclidean space with Sommerfeld's 
radiation condition we focus on the Helmholtz problem with Robin boundary condition on 
the bounded open set X. Problem (1.1) is well-posed in H X {X) for all g G L 2 (dX). In fact, 
writing a variational formulation of (1.1) shows the uniqueness of a solution to (1.1), while 
the existence of a solution follows from Fredholm's alternative. 

The thermo-acoustic imaging problem can be formulated as the inverse problem of 
reconstructing the absorption coefficient q from thermo-acoustic measurements q\u\ 2 in X. 
The quantity q\u\ 2 in X is the heat energy due to the absorption distribution q. It generates 
an acoustic wave propagating inside the medium. Finding the initial data in the acoustic 
wave from boundary measurements yields the heat energy distribution. Our aim in this 
paper is to separate q from u. We provide an explicit formula for reconstructing q from the 
heat energy q\u\ 2 in X. As far as we know, our formula is new. Indeed, it is promising since 
it can be used as an initial guess to achieve a resolved image of the absorption distribution 
in a robust way. 

Our first task is to enrich the set of data. Suppose that we have measurements q(x)\uj\ 2 
corresponding to linear combinations of boundary data gj, for j = 1, . . . , d + 1. We show 
that one can construct the set of quantities: 

£ = {Ej(x) = q(x)uj(x)ui(x), x G X | j = 1, . . . , d + 1}, (1.2) 
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where Uj denotes the solution of 

(A + k 2 + ikq)ui =0, x £ X, , x 

v (1.3) 
• Viij — ikuj = gj, x £ (9X, 

provided that (fl^ji* is a proper set of measurements (see Definition 2.1). The construction 
of Ei was completely described in [12] and that of Ej, j = 2, . . . , d + 1, will be done using 
Proposition 2.6. Noting that 

Ui Ei 

— = ~ET, j = 2,...,d+l, 

we are able to establish an exact formula for q provided that £ = ^ s "g°°d" enough 

as in Theorem 3.3. This procedure will be described in Section 3. 

As said, the collected data £ are often corrupted by measurement noise that varies on 
very small length scale. This renders the aforementioned exact formula, which requires 
differentiating the data up to third order, completely unpractical. To solve this issue, we 
smooth the noise by averaging the data over a small window and apply the smoothed data 
to the exact formula. The resulting function is then shown to be close to the real one, 
provided that the width of the averaging window is properly chosen. We thus view this 
function as an initial guess and then perform a further step of least square optimization. 
The resulting reconstruction improves the initial guess in the 1? sense. 

The rest of the paper is organized as follows. In Section 2 we introduce the notion of 
a proper set of measurements and its role to get data £ and some useful estimates as well. 
The aim of Section 3 is to provide an explicit formula for reconstructing q when a proper 
set of measurements is given. In Section 4 we study the Frechet differentiability of the 
data with respect to variations of q and prove that the differential operator is invertible for 
small enough variations. In Section 5 we consider a noise model for the data and show how 
to regularize the exact inversion formula in order to obtain a good initial guess. We also 
perform a refinement of the initial guess using an optimal control approach and show that 
this procedure yields a resolution enhancement. 



2 Preliminaries 

Motivated by [6], we introduce the following concept. 

Definition 2.1. The set (gj)jt.\ C L 2 (dX) is a proper set of measurements of (1.1) if 
and only if: 

(i) \u\\ > in X. 

(ii) The matrix [uj, V %]i<7<d+i is invertible for all x £ X. 

Here, T denotes the transpose and Uj is the solution of (1.3). 

The following proposition is a direct consequence of Lemma 4.1 in [12] and Proposition 
3.1 in [11]. It plays an important role to prove that it is possible to find a proper set of 
measurements. 
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Proposition 2.2. Let 5 > and m > d/2. There exists a positive constant C such that for 
any £ £ C rf ,£ ■ £ = 0, and |£| > (5, and for any q S H m (X), the solution w of 

&w + £-Vw = -{k 2 + ikq X {X)){l + w) in]R d , (2.1) 

where x(X) denotes the characteristic function of X, satisfies 

C\\q\\ H m(x) , . 

\\w\\h^(x) < j|j • (2.2) 

Proposition 2.3. If q £ H m (X), m> 1 + i/ien /tas a proper set of measurements. 

Proof. Let e be a small number. By choosing £ such that £ • £ = and |£| is large enough, 
we find from the Sobolev embedding theorem and (2.2) that the solution w of (2.1) satisfies 

IMU°°(X) + < e. (2.3) 

It is not hard to verify that the function 

u = e ix (l + w) 

is a solution of 

(A + k 2 + ikq X {X))u = (2.4) 

and it satisfies 



> le^l(l-e) > 0. 



Choosing g\ = v-Vu — iku on dX gives a solution u\ of (1.3) satisfying part (i) of Definition 
2.1. 



Define 



and 



£j = K e J +* e i+l)> J = 1, • . - , <i — 1, 
£d = n(ed + iei), 

✓ d-1 d-1 x 

£ d +i = n ( J" ej + Vd - le d + i \ ^ e i - Vd - le d J 

^ .7 = 1 .7=1 ' 



j=i j= 

where n>l and ej is the jth component of the natural basis of M. d . Again, it is not hard 
to verify that 

for all j = 1, ... , d+1, and the vectors (1, £j)i<j<d+i are linearly independent in C d . Hence, 

det [! ej] 1 < y < d+1 |»l, (2-5) 
provided that n>l, Let tOj , 1 < j < d + 1 , be the solution of 

Awj + 2£j ■ Vwj = -(k 2 + ikqx(X))(l + Wj) 
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and 

i 

be the solution of (2.4). We have 



det [ uj V T Uj ] ^ d+1 = e^ x (l + Wj ) det (1 + Wj ) £j + 



Thus, (2.3), (2.5), the continuity of the map that sends a square matrix to its determinant 
and the choice of large n imply the second part of Definition 2.1 with 

gj = v • Vuj — ikuj, j = 1, . . . , d, 

on dX. □ 

Remark 2.4. The solution w of (2.1) is the so-called complex geometric optics solution of 
(1.1), which was introduced in [10, 29]. The proof of Proposition 2.3 was partly motivated 
by [30]. 

We next construct the data £, mentioned in Section 1. Let us for the moment accept 
the following proposition. 

Proposition 2.5. If g is given, then one can make some measurements to obtain q(x)\u\ 2 , 
x £ X, where u solves (1-1). 

The following proposition holds. 

Proposition 2.6. Let gi,g2 £ L 2 (dX). Denote by Uj the solution of 

(A + k 2 + ikq)ui = 0, x e X, , x 

v j = 1,2. (2.6) 

v • Vuj — ikuj = gj, x E dX, 



Then the function q(x)u2(x)ui(x),x G X can be evaluated. 

Proof. Applying Proposition 2.5 for g\ + g2 and then ig\ + g2, we obtain the knowledge of 

q\u\ + U2\ 2 and q\iu\ + U2\ 2 , 
respectively. Then the desired data E2 is given by 

E 2 = ^(q\ n i +u 2 \ 2 - <?|ui| 2 - q\u 2 \ 2 ) + ^{q\iui + u 2 \ 2 - q\ui\ 2 - <?|ti 2 | 2 ), (2.7) 

which can be easily verified. □ 

Let (gj)1j=l be a proper set of measurements of (1.1) and Uj be the solution of (1.1) with 
g replaced by gj. From now on, we have the knowledge of 

£ = (Ej)<*±l (2.8) 

where Ej = qu{Uj, and £ is, therefore, considered as the data to reconstruct q. 

We also need the following proposition. It plays an important role to evaluate the 
derivative of the data with respect to q in Section 4 as well as some crucial properties. 
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Proposition 2.7. Let q £ L°°(X) be such that inf q > 0. For a// / £ L 2 (X), £/ie problem 

(A + k 2 +ikq)u = f, xeX, 
v ■ Vn — z/cu = 0, x £ <9X, 

/ios a unique solution. Moreover, the solution satisfies 

Ml2 ^ - khr q mL2 ^ (2 - 10) 

and 

V(fc2 + 1)+A;inf g |m| 

IMItf 1 ^) ^ ^-f^ Wfh^x)- (2.H) 

Proof. The well-posedness of (2.9) is well-known. Using the test function u in (2.9) and 
considering the imaginary and real parts of the resulting equation, we can establish (2.10) 
and (2.11), respectively. □ 



3 The exact formula 

The main aim of this section is to reconstruct q when a proper set of measurements (<7j)jiJ 
of (1.1) and the data 8, defined in (2.8), are given. 
Let 



«i = §-, 2<j<d + l. (3.1) 



Then it is not hard to see that 

Uj = OCjUl, 

for 2 < j < d + 1. We have the following lemma. 

Lemma 3.1. Let j3 = 9(uiVui). Then 

-div/3 = kE u in X. (3.2) 

Proof. Let ip £ C^°(X, M) be an arbitrary function. Then using (pu\ £ Hq(X) as a test 
function in 

— Au\ = (k 2 + ikq)u\ 

yields 

/ <f\Vu\\ 2 dx+ / tii Vui ■ V(pdx = / (k 2 + ikq)\u\ 2 (fdx . 
Jx Jx Jx 

Taking the imaginary part of the equation above gives 

— / div (9uiV«i)(/?dx = / kq\u\\ 2 ifdx = / kE\<pdx, 
Jx Jx Jx 

and (3.2) follows. □ 

The following lemma plays an important role in the derivation of an exact inversion 
formula for q. 
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Lemma 3.2. For all 2 < j < d + 1, 

Va,.(viogJ--^M)=Aa J , (3.3) 

Proof. Let us fix j 6 {2, . . . ,d+ 1}. Since iij is a solution of the Helmholtz equation under 
consideration, 

(k 2 + ikq)ocjUi = — A(ajUi) 

= — otjAui — u\Aaj — 2Vui • Vay 

= (fc 2 + ikq)a.jU\ — u\Aotj — 2Vu\ • V«j. 

Therefore, 

— E\Aaj = 2quiVui ■ Vaj 

= q (V|ui| 2 + 2i$PujVu 1 ) ■ Voy. 

We have proved that 

-E 1 Aa j = q (V\ui\ 2 + 2i0) ■ Voy, 

or equivalently, 

qV\ui\ 2 ■ Vaj = -E 1 Aa j - 2iqf3 ■ Vaj. (3.4) 
On the other hand, differentiating the equation E\ = <?|ui| 2 gives 

VE 1 = qV\ui\ 2 + Erf log q. 

This, together with (3.4), implies 

(V-Ei - SiVlogg) • Vaj = —EiAaj - 2iq/3 ■ Vay, 

and (3.3), therefore, holds. □ 
We claim that the set 

is linearly independent for all x £ X, where otj was defined in (3.1). We only prove this 
when d = 2. The proof when d is larger than 2 can be done in the same manner. In fact, 
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the linear independence of {Va2, VCK3} comes from the following calculation: 
det 



V T a 2 
V T a 3 



1 , 

— det 



u 



1 



det 



U\V T U2 — U2V T U\ 
MlV T ti3 — U^V T Ul 

U\V T U2 
MlV T U3 — UsV T U\ 



— ^ ( u\ det 



V T u 2 



— ? det 



ui V T ^i 
u 2 V T w 2 
^3 V T -u 3 



— U2 det 
+ U3 det 
— U2 det 



V T u 2 

V T Ul 



V 1 u 3 



Here, part (ii) in Definition 2.1 has been used. Since the d x d matrix 

A = [V T aj + i]i<j< d , Aji = diatj+i, 
is invertible, we can solve system (3.3) to get 

vlog ^-^T = G ' 



(3.5) 



(3.6) 



where a is the vector a = A~ l [(V T A T ) T ]. 

We are now ready to evaluate q. We first split the real and the imaginary parts of (3.6) 
to get 

Vlog — = — -V log Si = R(a) 
Ei q 



and 



Then, differentiating (3.8), we have 



2q 



(3.7) 
(3.8) 



div/3 



Si9(a) ■ Vg _ div (Si3(a)) 
2q 2 Tq ' 

This, together with (3.2) and (3.7), implies 

El (8(a) + Vlog Si) • 9(a) - div (f?i3f(o)) 
q ~ 2kW\ 



+ 



Sijft(a)- 9(a) + VS X - 9 (a) , Si div 9(a) + VEi • 9(a 
2fcSi 

»(a) -9(a) -div 9(a) 



2£Si 



2fc 



The results above are summarized in the following theorem. 
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Theorem 3.3. Given a proper set of measurements (g^^i so that the matrix A, defined 
in (3.5), is known and invertible. Then, 

-K(a) • 9(a) + div$s(a) 
Q{x) = ^ , (3.9) 

where a = ^ _1 [(V T ^ T ) T ] and A = {dia j+ i)jj =lj ... !d . 

Remark 3.4. Although in the proof of Theorem 3.3, we wrote some notations requiring the 
first and second derivatives of £ at a single point x € X, it is not necessary to impose the 
smoothness conditions for E. The reason is that one can make the arguments and establish 
(3.3) in the weak sense. We argued, using strong forms of differential equations, only for 
simplicity. 

Remark 3.5. Formula (3.9) is unstable in the sense that if there are some noises occurring 
when we measure the data Ej, 1 < j < d + 1, then q, given by (3.9), might be far away 
from the actual q since the right-hand side of (3.9) depends on the derivatives of the noise 
(up to the third order). 

4 The differentiability of the data map and its inverse 

Let < q min < q max . Let 

L°?{X) = Le L°°(X) : q min <p< q m ^ in x) . 



Then, L+(X) is an open set in L°°(X). We define the solution and the data map as 

u:L°°(X) -> H\X) 
q i—)- u[q] 



and 

F : L°°{X) L?(X) 

q ^ F[q]=q\u[q]\ 



2 (4-2) 



where u[q] is the solution of (1.1). The map F is well-defined because of the Sobolev 
embedding theorems and the fact that d = 2 or 3, which guarantees that u G L 4 (X). 

The main purpose of this section is to study the differential operator, T>F[q], of F and 
show that it is invertible provided that q ma , x is small enough. 

Lemma 4.1. The map u, defined in (4-1), is Frechet differentiable in L°^(X). Rs derivative 
at the function q is given by 

Vu[q](p)=v(p), V P £B q , (4.3) 

where B q C L°°(X) is an open neighborhood of q in L°°(X) and v{p) is the solution of 

(A + k 2 + ikq)v = —ikpu[q], ^ 
v ■ Vv — ikv = 0, x S dX. 
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Consequently, F is also Frechet differentiable and 

VF[q]p = p\u[q}\ 2 + 2q^(u[q}v(p)), Vg G L?(X), p G B q . (4.5) 



Proof. It is sufficient to show that 



lim h(p) = 0, (4.6) 
IIpIIl°°(x)— >o 

where 

s IHq + p] -«[?] - v{p)\\ L 2 {X ) 
Hp) = o • 

\\P\\l°°{x) 

In fact, since u[q + p] — u[q] — v(p) solves the problem 

(A + k 2 + ikq)(u[q + p] — u[q] — v(p)) = —ikp(u[q + p] — u[q]), x G X, 
v ■ V(u[q + p] - u[q] - v(p)) - ik(u[q + p] - u[q] - v(p)) =0, x G dX, 

we can apply inequality (2.10) to obtain 

II r j. 1 II < Ml <r \\p\\L~(.X)\\{u\q + p\ - u[q])\\ L 2 (x) 

\\u[q + p\- u[q\ - v{p)\\ L 2 {x) < ^— . (4.7) 

On the other hand, since u[q + p] — u[q] satisfies 

(A + k 2 + ik(q + p))(u[q + p] - u[q\) = -ikpu[q], x G X, 
v ■ (Vu[q + p] — u[q]) - ik(u[q + p] — u[q]) =0, x G dX, 

inequality (2.10), again, implies 

\\u[q + p \-u[q L 2 (x) < . . 4.8) 

v ; mf(g + / o) 

Combining (4.7) and (4.8) yields (4.6). Using the chain rule in differentiation, we readily 
get (4.5). □ 

Using regularity theory, we see that u[q] belongs to L°°(X) in the two-dimensional 
case. In three dimensions, we should assume that g G H 1 l 2 {dX) in order to claim that 
u[q] G L°°(X). Hence, T>F[q] can be extended so that its domain is L 2 (X). By abuse of 
notation, we denote the extended operator still by T>F[q\. The following key lemma of this 
section establishes an estimate of the L 2 (X) norm of v (/)), the solution to (4.4), in terms 
of the L 2 (X) norm of the source pu[q]. A corollary of this result allows us to show the 
invertibility of VF[q] from L 2 (X) to L 2 (X). 

Lemma 4.2. Assume that the origin is included in X and define 

iad(X) = sup \x\. 

x&dX 
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Suppose that X is star-shaped and balanced with respect to the origin so that 

x ■ v x > 7rad(X) 

for some positive number 7. If 



and k > 2, then 
where 



\q\\ L °oi&d(X) < i 



\v(p)\\ L 2 < r]\\pu[q]\\ L 2, 



i] 



(1 + 7- 1 ) 2 + 2d + 29 



' (11 - 2d) 

Proof. Let us define the bilinear form 

B[v, w] = 

and the linear form 



max{rad(X), 1}. 



A' 



Vv ■ Vwdx + k 2 I vwdx + ik / qvwdx + ik / vwds, 

dX 



X 



X 



G[w] = — ikpu[q]wdx. 
Jx 

Then the weak solution of (4.4) is characterized by v satisfying 

B[v,w] = G[w], ^weH 1 (X). 



(4.9) 
(4.10) 

(4.11) 
(4.12) 

(4.13) 



Using w = v in (4.13) and considering the imaginary and real parts separately, we have 



t)X 



x 



/ q\v\ dx < 


J 


' puvdx 




Jx 




X 




2 / \v \ 2 dx < k 


/ puvd 


r 


Jx 




Jx 





(4.14) 



It follows from these inequalities that 



~ L i(dX) ^ \\P u \\l 2 \\v\\l 2 , 



and 



k 2 

\\Vv\\h <(k 2 + l)\\v\\ 2 L2 + T \\ pu \\ 2 L2 . 
To estimate |M|l2, we mimic the technique used in [24, Chapter 8]. We have 

R(Vu ■ V(a? ■ Vv)) = \Vv\ 2 + x ■ v(^-), »(«(x-Vu)) =x- v(^j-). 
Integrating the first equation above gives 

/ R(VvV(x -Vv)) dx = [ \Vv\ 2 dx+- [ V • {x\Vv\ 2 ) - (V • x)\Vv\ 2 dx 
Jx Jx 2 J x 

= \f {v . x )\W\ 2 d S + (l- d -)\\Vv\\ 2 L2 . 

1 JdX 1 



(4.15) 
(4.16) 
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The second term above is due to the fact that V • x = d. Similarly, 

k 2 [ $t(v(x • Vu)) dx = ^ / V- (x\v\ 2 ) - {V-x)\v\ 2 dx 
Jx 2 J x 

= — — — f \v\ 2 dx H f (y ■ x)\v \ 2 ds. 

2 Jx 2 J ax 



Consequently, taking w = —x ■ Vf in (4.11) we find 

,2 -i r 

(x ■ v)\v \ 2 ds 



rile 2 1 f h l r 

-$tB[v,x-Vv] = —\\v\\ 2 L2 + - / {x-u)\Vv\ 2 ds- — \ 

+(! - ^)II Vu IIl2 + K ( -ifc [ qv(x-Vv)dx-ik f v(x-Vv)dt 
2 \ Jx Jax 

Equate the above expression with the real part of — KG[x • Vt>], i.e., Kifc J pux ■ Vv dx. We 
then obtain the estimate (using the fact that x ■ v > 7rad(X)): 

dfc 2 Il9 rad(X)7,,„ „, £; 2 rad(X) .. ,, 9 ,d . „„ ll9 

— Ikllia + v 2 ; i v^i| 2(ax) < — ^Mmi! w + (g - i)liv«ii| 2 

+/c rad(X) (|M|l°°|MIl2||Vv||l2 + ||«|Ua(ex) II ^Ah*{dX) + ||pu|| £ a||Vt;|| £ a) . 
On the other hand, it follows from Young's inequality that 

IMlL2(ex)l|Vv|| L 2 (ax) < e||Vt>|| 2 2(ax) + j- e \\v\\ 2 L 2 {gx) , 

for all e > 0. We choose e such that ke = 7/2 to get 



/c 2 rad(X) .. 

7rad(X) 2 k 2 iad(X) 7 + 1 



2 



Mll^ax) + fc r ad(X)||t;|| L 2 (ox) ||Vz;|| L 2 (9x) 



< 



■W Vv hHdx) + « — IMIl^x)- 



2 11 2 7 

Recall (4.15). The left-hand side of the inequality above can be further bounded by 

2 ^|v.li, ( «, + ^^2±i(«l|.|li, + ^IWIi.). ("7) 

Applying Young's inequality to the term ||/9w||l2 ||Vu||x,2 with efcrad(X) = 1/8 yields 

kmd(X)\\pu\\ L2 \\Vv\\ L2 < |||V^||| 2 + 2A: 2 rad 2 (^)|| /0 ^||| 2 . (4.18) 

Applying the same technique to the term ||v||x,2 ||V«||i2 shows 

krad(X)\\q\\ LOO \\v\\v\\Vv\\v < ^||V«||| 2 + 2A: 2 rad 2 (X)|| (? ||| 0O |M| 2 2 . (4.19) 

Finally, recalling estimate (4.16) and combining the above inequalities, we have 
d 2 /rad(X) 7 + l d 3 



7 ei + (| - |)(1 + k~ 2 ) + 2|| (? || 2 co rad 2 (X)) \\v\\l 2 

+ (£^l^ + 2rad 2 ( x ) + i ( ^^)| W l!2. 
12 



Suppose that the wave number k is larger than 2 and the product ||g||i<x>rad(X) is 
smaller than 1/4. Then, if 4ei is chosen to be (rad(X)(7 + , the coefficient in front 

of II v\\^2 on the right is less than 5d/8 — 11/16. Then ||u||£2 term on the left dominates and 
we have 

" - 8 )Mh < (^^rad 2 (X) + 2rad 2 (X) + \(± - ||p.||| 2 . 
Estimate (4.9) follows from this immediately. □ 

Lemma 4.3. Let -q denote the constant (4.10). Suppose that the absorption coefficient q is 
such that 

V\W\\l°°{x) < 1- (4-21) 

Suppose also that \u\q\\ is bounded from below by a positive number. Then the map T>F[q], 
as an operator from L 2 (X) to L 2 (X), is invertible. Moreover, 

\\DF[ar l \\ c{L 2 {x)) < 1 (4.22) 

mf \u[q]\ 2 Jl -4n\\q\\ L «>ix) 



Proof. Define 

T[q}(p) = \u[q}r 2 VF[q](p)-p. 

It is not hard to see that T is compact since it can be decomposed as 

T:L 2 (X) -> H\X) L 2 {X) -> L 2 {X) 

p i — y v(p) i — y v(p) i->- 2q\u[q]\~ 2 R(u[q]v(p)). 

The continuity of maps in the diagram above can be deduced from Proposition 2.7 and the 
choice of g such that \u[q]\ > in X. 

On the other hand, a straightforward calculation shows that 

\\VF[q}(p)\\ 2 LHx) > inf |n[(?]| 4 ||p||| 2(x) [l - 4r ? ||g||^ (x) ] . (4.23) 

In fact, 

\pF[l\(p)\\h(x) = I [p 2 |«[g]| 4 + W(«[g]lJ(p))+4gp|ub]| 2 ^(u[#(p))]^ 



x 



> influf 



x 



4q 2 R 2 (u[q]v(p)) 
u[q}\< 



pHlW + i., r „ii2 + *q&(pu[q]v(p)) 



dx 



> mi\u[q]\ 2 J [p 2 \u[q]\ 2 - 4\\q\\ Lx(X )\pu[q]v(p)\] dx 

> inf \u[q]\ 2 \\p\u[q]\\\ 2 L2{x) - 4||g|| Z) cx> (X ) \\pu[q] \\ L 2 {X ) h(p)\\tf{X) 

> mf\u[q]\ 2 \\pu[q]\\ 2 L2{x) [l - 4r]\\q\\ Lxi{x) ] . 

Since ?7||g||£oopQ < 1/4, we find (4.23). It follows that the kernel of T>F[q] is {0}. Hence, 
by the Fredholm theory, T>F[q] is invertible. Moreover, (4.23) also implies (4.22). □ 

Remark 4.4. Recall the definition of 7] in (4.10). When X is a ball, r\ is roughly three to 
four times the radius of X in dimensions three or two. Condition (4.21) hence requires that 
||g||i,oorad(X), which can be interpreted as the typical absorption rate as signals propagate 
to the boundary, should be sufficiently small. 
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5 Measurement noise and resolution enhancement 



In this section, we consider additive noise in the data set £ given in (1.2). 



5.1 Noise model 

As described in Proposition 2.6, the data £ are acquired by measuring several sets of ab- 
sorbed radiations: q\u\ + Uj\ 2 , q\iu\ + Uj\ 2 , q\u\\ 2 , and </|ttj| 2 for j = 2, . . . , d+1. In practice, 
the measurements of these absorbed energies are corrupted by additive noises. We model a 
typical energy measurement by 

E m {x) = E(x) + aW s (x). (5.1) 

Here and in the sequel, the superscript "m" indicates measured quantity, and E itself is the 
pure quantity without noise. W$ is a stationary random field with mean zero and covariance 
function of the form 

E[W s (x)W s (y)] =E[Ws(0)Ws(x-y)] =i?(^Q, (5.2) 

where R is an integrable function normalized so that R(0) = 1. In this additive noise model, 
a 2 is the variance of the noise, 5 is the correlation length which is related to the distance 
between measurement points. 

The random process W$ is assumed to be bounded almost surely by a constant inde- 
pendent of 5. This constant is assumed to be smaller than E m i n which is a lower bound 
for the real energy. This technical hypothesis ensures that E m is bounded from below by a 
positive constant for any a < 1 and for any 5. 

In the forthcoming analysis, both the noise variance a and the noise correlation length 5 
will be supposed to be small. We assume that the measured data £ m = (EJ 1 )^^ are given 
by 

E?{x) = E 1 (x) + aW 51 (x), 

Ef{x) = Ej(x) + aU Sj (x) + iaVsj(x), j = 2, . . . , d + 1. 

According to the procedure of measuring Ej, the random fields Ugj and V$j are given 
by (Wsij — Wsi - Wgj)/2 and {Wgij' — W§i — W$j)/2 respectively, where Wgj, Wgij and 
Wsy correspond to the additive noises of the energy measurements g|itj| 2 , q\u\ + Uj\ 2 and 
q\iu\+Uj\ 2 , respectively. It is natural to assume that W$i, W$j, Wgij and Wsiy are mutually 
independent and have the same statistical distribution as Wg in (5.1). As a consequence, 
U$j, V$j and W$\ are correlated. 



5.2 Initial guess with smoothed data 

We smooth the data £ by using the convolution kernel 

where p € (0, ^g) and ip is in the Schwartz space of smooth nonnegative functions that 
decay rapidly at infinity and that satisfy f^ d (p(x)dx = 1. The condition p < d/(d + 6) will 
be clear later. The following lemma will be useful. 
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Lemma 5.1. Let \ j\ denote the sum of all components of the multi-index 7 and d^ip (resp. 
d^ips) denotes the usual 7— partial derivative of if (resp. (ps). For any 5 we have 

E|W w *aV 5 | 2 <^ (d+2|7l)p ||i2||i 1(Rd )||5V||| 2(Rti) . (5.5) 
More precisely, for i5<1, we have 

E\W ls *d~ f vs\ 2 = S d - {d+2hl)p [ R{y)dy [ \d^p{y')\ 2 dy' + (5 d ' {d+2hl)p ) . (5.6) 

JR d jR d 

Proof. The variance (5.5) can be written as 

E\W ls *<P<p s \ 2 =E * +2dp / [ W ls (x-y)W ls (x-y')(d^)(^)(d^)(^)dydy' 

R{ y^ ){d ^ ){ y ){d i^ )dydy '. 



S2 P \l\+2dp J Rd J Rd V J A tVVjpM 

We apply the change of variable (y — y') / 5 — )• y' and y/<P — )• y, and take advantage of the 
resulting Jacobian. We verify that the variance can be written as 

E\W 15 *d^s\ 2 = S d+dp - 2phl - 2dp [ R(y') [ (5»(y)(0»(y - S^y^dydy' . 

JR d JR d 

Using Cauchy-Schwarz inequality and the fact that d^tp G L 2 and R G L 1 , we obtain (5.5). 
Since cTV £ -^ 2 > P < 1> an d is integrable, (5.6) is also easily verified by the dominated 
convergence theorem. □ 

Remark 5.2. The above calculation works also for Ujs and Vjg. 

We smooth the data by evaluating the convolution with the kernel (ps- 

E^Ef*^ j = l,...,d + l, (5.7) 

which gives 

E\ = E 1 *ips + o-Wis*<p S , (5.8) 

Ej = E j *ips + o-U jS *ips + io-V j s*(p6, j = 2, ...,d + l. (5.9) 

Here and below, the superscript "s" indicates smoothed quantities. The parameter 5 P can 
be interpreted as the size of the averaging window. To simplify the notation, Ej§ will be 
used as the short-hand notation for the smoothed unperturbed data Ej * (ps in the sequel. 

Proposition 5.3. If we substitute the smoothed measured data (E*) d ~^\ into the reconstruc- 
tion formula (3.9); 

,-(,) = + *■*(■■•), (5 . 10) 

with a s = (A s y l [(V T (A s ) T ) T ], A s = (dio^ j+1 ) jt i =lt ... 4 , and a) = EpEf, then the estimate 
q s satisfies: 

supE[|g s (x) -q S (x)\ 2 ] < Ca 2 5 d ^ d+ ^ p , (5.11) 
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where 

-$t(a 5 )%(a s ) + div%{a s 



2k 

\d+l 



is obtained by substituting the smoothed unperturbed data (Ejs)j^ 1 into the reconstruction 
formula (3.9). 

Proof. We substitute the smoothed data (-Ey)jiJ into the reconstruction formula (3.9). 
Recall the definitions of A and ay in (3.5) and (3.1). Then, 



a s = E jS + (rUjS *<P6 + ivVjS * <P6 (5 u) 



When a <C 1, we can linearize this term and find that 

a] = %i - + + io y^™ + {a 2 ). (5.13) 

Eu Eu E 1S Eu Eu 

The coefficients of the matrix A s are defined by A s ^ = dia S j +1 and they can be expanded 
from (5.12) as 

A s fl = A 5 jt + aAf ] + o(a5 d ^ d+2 ^ 2 ), l<j,l<d, (5.14) 

where 

A s _ q E j+is A 5(i) _ Wu* djips Ej+u | Uj+is*diips | yj+is*dnps 
3 Eu ' jl E 1S Eu E 1S Eu 

The leading-order error terms o~A S ^"' have zero means and their variances are of order 
0(a 2 5 d ~( d+2 ^ p ) according to Lemma 5.1, provided that the functions Ej's are sufficiently 
smooth with bounded derivatives. The following error terms like Wu * <-Psdi ( E ti is ) are 

smaller since their square means are of order 0(a 2 5 d ~ dp ). 

bmce A 5 is a smoothed version of A, which was defined in (3.5) and whose determinant 
can be bounded from below by a large constant (see Proposition 2.3), the inverse of A s is 
well defined. Linearizing (A s ) _1 , we have 

(^r 1 = {A 5 )- 1 +a{A s r 1 A^ 1 \A s r 1 + o(a5 d / 2 ^ d+2 ^ 2 ). 

Similarly, the vector (V T A sT ) T can be decomposed as 

cy7T,sTn _ (T7 t a t\ ,„( Wu*Aip s E j+ u Uj+u*A^s , Vj+u*A(ps 

l v A )j ~l v A + p p 1 p ri 

V Eu Eu Eu Eu 

+ o(a5 d ' 2 ^ d+A ^ 2 ). 
Finally, we have for the vector a s = (A**^ 1 (X7 T A sT ) T : 

s , STfA8\-if Wis * Atps Ei+u . Ui +ls *A(p 5 ,%j*A w 

a,- =a„- + a > (A ) „■> — — — \- i — 

3 3 ^ 3 V E l6 Eu Eu Eu 

+ o(a5 d / 2 -( d+4 ^ 2 ), 
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and 

diva* =diva 5 + a V (A 6 )', 1 f.MMi^jj + + ^+1**0^)^ 

3,1=1 

+ o(a5 d / 2 -^ p / 2 ). 

The vector a* 5 = (A s ) 1 (\7 T A sT ) T is obtained by applying formulas (3.1) and (3.5) to the 
smoothed unperturbed data (Ejg) d ~^. The leading-order error terms have zero means and 
their variances are of order 0(a 2 5 d ~( d+ ^ p ) according to Lemma 5.1. Our choice p < 
guarantees that the noisy data are smoothed enough so that the terms above have variance 
of order smaller than a 2 . To summarize, if we apply (3.9) to the smoothed data (Ej)-t^, 
then we get 

fix) -ZrlY S(A>tf (_ w ^*d^ 5 E^s + U l+15 * dj A<p s 



2k\ ^ v Jjl V Eis Eu E u 

K 3\i=i V (5.15) 

+ M(A s )- l Vl+1S * djAlps ) + o(aS d / 2 ^ d+ ^/ 2 ), 



E 15 

from which we deduce the desired result. □ 

The terms q$ can be shown to be close to the real absorption parameter q Q uniformly in 
x (we show this in Theorem 5.4). However, it is impossible to separate qs from the noise, 
that is the other terms in (5.15). Nevertheless, the estimate q s is a good initial guess in the 
mean square sense as shown by the following theorem. 

Theorem 5.4. Suppose that the pure data (Ej) d +l belong to C 3 ' e for some positive real 
number e. Then, we have 

\\qs-qo\\ L °°(X)<C5 £ P. (5.16) 
As a result, estimate (5.10) obtained from the smoothed data satisfies 

supE[|g s (x) -q {x)\ 2 } < C(5 2£p + o- 2 5 d ~ {d+6)p ). (5.17) 

Proof. Under the conditions of the theorem, the inequality \SPEj{x — y) — d" 1 Ej{x)\ < C\y\ £ 
holds for some constant C and for any multi-index 7 with I7I < 3. As a result, we have the 
following estimate as an analog of Lemma 5.1: 



\d^E jS {x) -d^Ej(x) 



)dy 



< c jr P I d \y\ £ W{^)\dy = cs^ 



(5.18) 



Then the estimate of qs follows because the reconstruction formula in (3.9) depends con- 
tinuously on the data and their derivatives. For the second estimate, we apply the triangle 
inequality and use the control of the stochastic terms in the linearization procedure. □ 

Remark 5.5. Estimate (5.17) is a bit over pessimistic. Indeed, it does not imply that q s is 
positive, which is a physical constraint for the absorption parameter. We will exploit this 
remark in the next section. 
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5.3 The optimization step and resolution enhancement 

Now we refine the above initial guess q s by an optimal control approach. We seek for the 
least square estimate of the discrepancy functional 

J[q} = ( \F[q](x) - E\{x)\ 2 dx. (5.19) 
Jx 

Here, E\ is the smoothed data (E\ + aWu) * tp$ and F[q] = q\ui[q]\ 2 is the absorbed heat 
energy with boundary condition g±. 

In Theorem 5.4, the initial guess q s is shown to be close to the true absorption coefficient. 
This allows us to approximate the integrand in the definition of J by its linearization around 
q s ; that is, 

J[q] « [ \VF[q s ](q - q s ) - b s \ 2 dx, (5.20) 
Jx 

where 6 s = Ef — F[q s ] is the residue. In the case when VF[q s ] is invertible from L 2 to L 2 , 
the least square solution of the approximate discrepancy functional is given by 

q, = q * + (VFtf])- 1 ?. (5.21) 

The following result shows that q* is a refinement of q s in the mean square sense (compared 
to Theorem 5.4). 

Theorem 5.6. Recall that q Q denotes the true absorption coefficient and assume that the 
condition in Theorem 5.4 holds. We have 

E[|k* - qo\\h(x)] = °^ 2£V + o- 2 5 d -( d+ V p ) + 0{5 2p + a 2 5 d ^). (5.22) 

Proof. From the definition of 6 s and E\, the residue can be expanded as 

b s = E x - F[q s ] + {E 1S - Ek) + aW ls * <p 5 . 

Since E\ = F[q ], the difference F[q Q ] — F[q s ] can be linearized as VF[q s ](q — q s )+o(q — q s ). 
This, together with (5.21), implies 

q* - q Q = {VF[q s \Y l {aW 15 * <p s + (E 15 - E x ) + o(q - q s )}. (5.23) 

Lemma 5.1 shows that aW±s * <fi$ has mean square of order cr 2 ^ 1- ^; the calculation in 
(5.18) shows that £^15 — E\ can be bounded uniformly by CS P ; the term q a — q s is also 
controlled in (5.17). Consequently, since T>F[q s ] has bounded inverse (see Lemma 4.3), the 
desired estimate holds. □ 

Remark 5.7. Assume that q Q is bounded from below and above by two known positive 
numbers § m i n and § max . Let 

= min{max{^,g min },g max } G [g min , g max ]. 

We can see that 

ll<?* - Qo\\l 2 (X) < Ik* - Qo\\l 2 (X)- 
We note that there is no guarantee that q* is positive, but the modified version g* is. In 
addition to this advantage, the estimate above shows that q* is a better approximation of 
q Q in comparison with q*. Further, the range of allows us to make iterations for further 
corrections. 
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Remark 5.8. Finally, we note that the above result also shows that the optimization step 
enhances the resolution. In fact, from (5.21) it follows that q* contains higher oscillations 
than q s and therefore, yields a more resolved approximation of q . 

6 Conclusion 

In this paper we have derived an exact reconstruction formula for the absorption coefficient 
from thermo-acoustic data associated with a proper set of measurements. Using a noise 
model for the data, we have regularized this formula in order to obtain a good initial 
guess. We have also performed a refinement of the initial guess using an optimal control 
approach and shown that this procedure reduces the occurring errors and yields a resolution 
enhancement. A challenging problem is to estimate analytically the resolution. It would 
be also very interesting to study the reconstruction problem in the case of incomplete 
measurements, where the thermal energy is known only on an open subset of the domain. 
The numerical implementation of our approach in this paper is the subject of forthcoming 
work, which will be published elsewhere. 
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